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Abstract 

When 0 < a < 1, the approximation for the Caputo derivative 

n 

y (a) (*) = r( 2 - a )h a Z a t )y ( x ~ kh ) + ° i h2 ~ a ) > 

where = l,a^ = (n — l) 1_a — n l ~ a and 

a[ a) =(k- l) 1 "" - 2k 1 ~ a + (k + I) 1 ” 0 , (k = 1 • • • , n - 1), 

has accuracy O (h 2 ~“). We use the expansion of Yfc=o k ° determine 
an approximation for the fractional integral of order 2 — a and the 
second order approximation for the Caputo derivative 

1 n 

y {a) ( x ) = r(2 _ a)h a Z - kh ) + 0 ( h2 ) - 

where <5^ = a^ for 2 < k < n, 

S o a> = cr o“ ) _ C(«- l)Yl a) = cr i“ ) + 2 C(« - 1),4 Q) = 4 a) -C(«- 1), 

and C( s ) is the Riemann zeta function. The numerical solutions of the 
fractional relaxation and subdiffusion equations are computed. 
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1 Introduction 


Fractional differential equations are used for modeling complex diffusion pro¬ 
cesses in science and engineering [1-5]. The Caputo fractional derivatives are 
important as a tool for describing nature as well as for their relation to inte¬ 
ger order derivatives and special functions. The Caputo derivative of order 
a, when 0 < a < 1, is defined as the convolution of the power function x~ a 
and the first derivative of the function on the interval [0, x] 


y^ a \x) = D a y(x ) 


d a y(x ) 
dx a 


1 

T(1 — a) 



y\q) 

(x-0 a 


da. 


When the function y(x) is defined on the interval (—oo,x], the lower limit 
of the integral in the definition of Caputo derivative is — oo. The Caputo 
derivative of the constant function 1 is zero, and 


D a D 1 - a y(x)=y\x). 


While the integer order derivatives describe the local behavior of a function, 
the fractional derivative y^ a \x) depends on the values of the function on 
the interval [0, x]. One approach for discretizing the Caputo derivative is 
to divide the interval to subintervals of small length and approximate the 
values of the function on each subinterval with a Lagrange polynomial. Let 
x n = nh and y n = y(x n ) = y(nh), where h > 0 is a small number. The 
Lagrange polynomial for the function y'(x) at the midpoint Xk- 0.5 of the 
interval [xk-i,Xk] is the value of y'(xk- 0 . 5 )■ 

Approximation ((TJ) for the Caputo fractional derivative is a commonly 
used approximation for numerical solutions of ordinary and partial fractional 
differential equations [6-8]. 


T(l-a)y^\x n ) = 


0 \ x n sj fc=1 J x k 

y(xk) - y(xk-i) f kh 


Xk 


y'{xk- 


0.5 j 


E 

k =1 

n 

E 

k =1 


1 


da 


'(fc-i)^ ( nh - ay 


:da 


y k - y k -1 ((n -k + 1 )hy 01 - ((n - k)h) 
h 1 — a 


\ 1 —a 
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Let p^ = (n — k + l) 1 Q — (n — k) 1 

n n n 

r (2 - ot)h a yi a) « ^p^iVk - S/fc-i) = Erf “ E^- 1 ^ 

fc=l fc=l /c=l 

71—1 

(a) , / (a) (a) \ (a) 

= Pn ’Vn + 2^ Vk [Pk - Pk-1) ~ Pi yo- 


Then 


k =1 


72 — 1 


O SS*) 
y n 


r(2 - a)h° 


Pn ] yn + E Vn-k ( Pn-k ~ Pn-k+l) ~ Pl^V 0 I ■ 


k =1 


Let ag Q) = pi“^ = 1, <Tn“ ) = — = (n — l) 1 a — n 1 a and 

d“ ) = - Pn \ + 1 = (* +1)‘-“ - 2fc‘-“ + (* - i) 1 -”, 

for k — 1, 2, - ■ • , n » 1. Denote 


r («) 


,(«) 


•A-hUn ^ ^ Vn—k- 
k=0 

We obtain the approximation for the Caputo derivative 

^ “ r(2-a)> 1 ° Aj '"- 

Approximation (pQ) has accuracy 0(h 2 ~ a ) when y 6 C^fO, x n ] (0). 


a) 


Table 1: Error and order of approximation (CQ) for y(a;) = cosx on the interval 
[0,1], when a = 0.6. 


h 

Error 

Ratio 

Order 

0.05 

0.0023484 

2.69618 

1.43092 

0.025 

0.000878437 

2.67338 

1.41867 

0.0125 

0.000330265 

2.65979 

1.41131 

0.00625 

0.000124548 

2.65171 

1.40692 

0.003125 

0.0000470549 

2.64687 

1.40429 
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The numbers a^ have the following properties: 

OO 

4 Q) > 0 , < 4 a) < ■■■ < a ( k a) < ■■■ < 0 , ^a ( k a) =0. 

k =0 

Approximation (0Q) and its modifications have been successfully used for 
numerical solutions of fractional differential equations, as well as in proofs 
of the convergence of numerical methods. One disadvantage of ([1]) is that 
when the order of the Caputo fractional derivative a ~ 1, its accuracy de¬ 
creases to 0(h). The numerical solutions of multidimensional partial frac¬ 
tional differential equations require a large number of computations, when 
the approximation has accuracy 0 (h). 

In section 4, we determine the second order approximation (J4|) for the 
Caputo derivative by modifying the first three coefficients of (jT]) with values 
of the Riemann zeta function. Approximation (HJ has accuracy O (h 2 ) for all 
values of a between 0 and 1. 

The ordinary fractional differential equation 

y W + By = F(t), (2) 


is called relaxation equation when 0 < a < 1, and oscillation equation when 
1 < a < 2. In section 5 we compare the numerical solutions for the relax¬ 
ation and the time-fractional subdiffusion equations for discretizations ([T]) 
and (J4|). We observe a noticeable improvement of the accuracy of the nu¬ 
merical solutions using approximation (J4]) for Caputo derivative, especially 
when a ~ 1. 

When y(x) is a sufficiently differentiable function, the integral in the 


definition of the Caputo derivative has a singularity at the point x. Sidi 10 
discusses approximations for integrals with singularities. 

The sum of the powers of the first n — 1 integers has expansion 11 


n—1 


n 


a+l 


^fc« = c (-«) + W t E(“ + 1 )% 

k= 1 m =0 v 7 


(3) 


where a ^ —1 and B m are the Bernoulli numbers. In section 3, we use 
expansion (J3]) to determine a second order approximation (|6|) for the left 
Riemann sums and the fractional integral of order 2 — a. In section 4 we 
determine the second order approximation for the Caputo derivative (|4|) from 
(J6]) , using discrete integration by parts and second order backward difference 
approximation for the second derivative. 
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2 Preliminaries 


In this section we introduce the basic definitions and facts used in the paper. 
The fractional integral of order a is defined as the convolution of the function 
y(x) and the power function x a ~ l on the interval [0,x] 



where a > 0. The fractional integral of order a is often denoted as y ( a \x ). 


The value of the fractional integral of order a of the constant function 1 is 
x a /T(a + 1). The Caputo derivative is defined as the composition of y'(x ) 
with a fractional integral of of order 1 — a. In Claim 1, we represent the 
Caputo derivative with the composition of the second derivative y"(x) and a 
fractional integral of order 2 — a 



The composition of fractional integrals satisfies 


J a jPy(x) = jPj a y(x) = J a+P y(x). 


The composition of the Caputo derivative and the fractional integral of order 


a, when 0 < a < 1, has properties 

D a J a y(x ) = y(x), J a D a y(x ) = y(x) - y( 0). 

In Theorem 3 we use the expansion of the sum of the powers of the first n — 1 


integers (|3jl . to determine the second order approximation for the fractional 
integral of order 2 — a 



where h = x/n, and £(s) is the Riemann zeta function, defined as the analytic 
continuation of the function 



n= 1 
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In the special case of (15|) . when a = — 1, the sums of the harmonic series 
have expansion Q 

n ~ 1 i 1 00 R 1 

E 1 , 1 \ ^ n 2m i 

k 2n ^ 2 m n 2m 

k =1 s=l 

where 7 ~ 0.5772 is the Euler-Mascheroni constant and i? 2 m are the Bernoulli 
numbers. 

I 11 section 4 we use the approximation for the fractional integral (J 6 lh to 
determine the second order approximations for the Caputo fractional deriva¬ 
tive 

yi ° >(x) = r(2 -a)h° g ^ ~ kh) - W^) y " {x)h ^° + ° b 2 ) ’ 
^^ = r(2 - a)h° E 5 ^ y ^ x ~ kW> + 0 M ■ ( 4 ) 

The numbers 5^ are computed from the coefficients of ([I]) by 

Si, a) = &o a) - CO - 1), <^i Q) = v{ a) + 2C(« - 1), 4 a) = 4 a) - CO - 1), 

4 a) = 4 Q) ( k = 2 , 3, - - - ,n). 

The values of the Riemann zeta function satisfy jl3j | 


C(*) = 


EhTTlE^ 1 ) 


1 _ 2 l ~ s ^ 2 n+1 

n =0 k =0 


k ‘ n k )(k + l)- s , 


for all s G C, and the functional equation 

CM = 2“7r“- 1 sin (^) r(l - s)C(l - s). 

From the functional equation for the Riemann zeta function we obtain a 
representation of ((a — 1)/T(2 — a) 


CQ-i) 

T(2 - a) 


= -2 


“-V^cos (™)c(2-a). 
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3 Approximation for the Fractional Integral of Order 2 — a 


In this section we determine a second order approximation (j 6 j) for the frac¬ 
tional integral of order 2 — a, when 0 < a < 1 

J 2 ~ a y(x) = —^-- [ {x - 0 1 _< WK- 

i (2 — a) J Q 

Approximation (JHD uses the left Riemann sums of a uniform partition of 
the interval [0, a:], and the values of y( 0) and y(x). The Caputo derivative 
y( a \x) = J 1 ~ a y'{x ) is dehned as the composition of the fractional integral 
of order 1 — a and the first derivative y'(x). In Claim 1 we use integration 
by parts to express the Caputo derivative as a composition of the fractional 
integral of order 2 — a and the second derivative y"(x). 

Claim 1. Let y G C 2 [0, x], and 0 < a < 1. 

T(2 - a)y {a \x ) = T(2 - a)J 2 ~ a y"(x) + y\0)x l ~ a . 


Proof. From the properties of the composition of fractional integrals and 
Caputo derivatives 

J 2 ~ a y"{x ) = (x) = J 1 ~ a (y'(x) - y'( 0)). 


Then 


y (a \x) = J^yfx) = J 2 ~ a y”(x ) + 



J 2 - a y"(x ) + 


y’{f))x l a 

T(2 — a) ' 


□ 


Let x = nh, where n is a positive integer. Consider the partition Vh of 
the interval [0, x\ to n subintervals of length h. Denote by C[^ h and the 
left Riemann sum and the Trapezoidal sum of the function (x — f) 1 ~ a y(f) 
for partition Vh 

n 7i—l 

d^l = h — mh) 1 ~ a y{mh ) = h nh — 

m =0 m =0 

/ ( nl 

Ty*h = | ( (n/i) 1 _a /( 0 ) + 2 ^(nh- mh) l ~ a y(mh) 

\ m=l 
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Substitute k = n — m 




k=1 

Xji = '^Y-x^h + h 2 ~ a Y kl ~ a y(x - kh ). 

k=1 

The numbers and Tyff are approximations for T(2 — a) X 2 ~°^y{x) and 

(5) 


r( a ) q-’i OL ) y(X) 1 —au 

S ,h 'yL ~ 2 


Now we use (J3]), to determine a second order approximation for the left 
Riemann sums of the constant function y(x) = 1. 

Lemma 2. Let x = nh, where n is a positive integer. 


4?j! = It- + + c(« - i)^ 2 " Q + O {h 2 ) ■ 


2—a j 


2 — a 2 
Proof. Consider the first terms of 


n —1 

k=1 

n 

E fc1 ' 

fc=i 

Multiply by /r 2_a 

n 

h 2 ~ a y kl ~ a 


n 2 a n 1 “ 


2 — a 


n 2 a n 1 Q 


+ C(« — 1) + ^ ( — ) j 


rr 


2 — a 


H-x-1- C(« — 1) + O — • 


rr 


x 


2—a 


k=1 


2 — a 2 


Ppx 1 a h + ((a — l)h 2 a + 0 


h 2 ~ 


n u 


We have that h 2 a /n a = h 2 /x a . Hence 

, % ry^—OL 1 

•»(o;) _ % | _ 1 — 


^-l~ px 1 a h + ((« - 1)L 2 “ + O (h 2 ) . 

2 — a 2 


□ 












In the next theorem we determine a second order approximation for the 
left Riemann sums of the fractional integral J 2 ~ a y(x ) when the function y(x) 
is a polynomial. 

Theorem 3. Let x — nh and y(x) be a polynomial. 

= r ( 2 - a ) j2 ~ Q y( x ) + y -^-x l - a h + C(« -1 )y(x)h 2 ~ a + O {h 2 ). (6) 

Proof. Let y(f) be a polynomial of degree m. The Taylor polynomial for y(f) 
of degree m at the point £ = x is equal to ?/(£). 

m 

V(0 =P0+Pl{x~f) + • •• +Pn{x -£) m =p 0 + ^2pk(x-f) k . 

k=\ 

Denote 

m 

s/o(0 = y(0 - y(x) = J2p k (x - £) fc . 

k =1 

The function {x — £) 1_c h/o(£) has a bounded derivative on the interval [0,x]. 
The trapezoidal approximation is a second order approximation for the 
fractional integral T(2 — o) J ( - 2 ^ a ' > yo(x). From Lemma 2 and (|5]) 

Pf = + C(<* - + O (IS) . 

2 — o 

Then 

Pt = Pt+PoPt = r( 2 -«) X-°y 0 (x)+ : X!L- x i-<‘+ P0 C(a-l)h 2 -“+O (A 2 ) . 

We have that y(x) = p 0 and J 2 ~ a 1 = x 2 ~ a /Y(?> — a), 

T(2 - a)J 2 ~ a y(x ) = T(2 - a) J 2 ~ a y 0 {x){x) + —x 2_ “. 

2 — a 

Hence 

7$ = = T(2 - a)J 2 -«y(x) + y(x)«a - 1 )A 2 "“ + O (A 2 ) . 

□ 
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In Theorem 3 we showed that (151) is a second order approximation for the 
left Ricmann sums and the fractional integral of order 2 —a, when the function 
y(x) is a polynomial. From the Weierstrass Approximation Theorem every 
sufficiently differentiable function and its derivatives on the interval [ 0 , x] are 
uniform limit of polynomials. The class of functions for which Theorem 3 
holds includes functions with bounded derivatives. In section 4, we present 
a proof for the second order approximation (J3J) of the Caputo derivative. 


Table 2: Error and order of approximation ((HD for y(x) = cosx (left) and 
y(x) = In (a; + 1) (right) on the interval [0,1], when a = 0.4. 


h 

Error 

Order 

h 

Error 

Order 

0.05 

0.00011853 

1.95822 

0.05 

0.00020451 

1.98580 

0.025 

0.00003019 

1.97331 

0.025 

0.00005145 

1.99083 

0.0125 

7.63 x 10 " 6 

1.98275 

0.0125 

0.00001292 

1.99400 

0.00625 

1.92 x 10 " 6 

1.98876 

0.00625 

3.24 x 10 " 6 

1.99606 

0.003125 

4.83 x 10 " 7 

1.99264 

0.003125 

8.11 x 10" 7 

1.99740 


4 Second Order Approximation for the Caputo Derivative 

In this section we use approximation (J 6 ]) to determine a second order dis¬ 
cretization for the Caputo derivative of order a, by modifying the first three 
coefficients of approximation ([!]) with the value of the Riemann zeta function 
at the point a — 1 . 

Denote by A \y n and A \y n the forward difference and the central differ¬ 
ence of the function y(x) at the point x n = nh. 

Vn- 1-1 Uni 

A h.yn Vn+ 1 2 y n y n — i- 

When y(x) is a sufficiently differentiable function 

<4+0.5 = % + O (V). y'; = ^ + o (h 2 ). 

Lemma 4. 

71—1 

A h y n = fcl ”“A \y n - k + n 1 _Q A^ 0 . 

k=l 
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Proof. 


n n—1 

AhVn = Y ^yn-k = ^Vn + X ^ Vn-k + ^Vo 
k= 0 k= 1 

n—1 

=Vn + Y(( k ~ 1 ) 1 ~“ “ + ( k + 1 ) 1 ”“) 2/n-fc + ^S/O 

fc = l 

n—1 n—1 n—1 

~Vn + X( fc - l) 1_a ?/„-fe - 2 X k 1 -°‘y n _k + Y( k + l f~ a Vn-k + ^ Q) yo- 

/c=l /c=l k=1 

Substitute K — k — 1 in the first sum and A' = k + 1 in the third sum 

n—2 n—1 n 

Ay n = Vn + Y Kl ~ a Vn-K -1 - 2 X kl ~ a Vn-k + X K 1 ~ a y n _ K+ 1 + (T^ a) y 0 . 

A'=l fe=l A'=2 

We have that 

n—2 n—1 

XI kl ~ a Vn-k-l = Y kl ~ a Vn-k- 1 - (n - l) 1_Q yO, 

fc=l /c=l 

n n n—1 

Vn H - ^ ^ k Vn—k+l ^ ^ ^ Vn—k+l ^ ^ ^ Vn—k-\-l y\. 

k=2 /c=l /c=l 

Then 

n—1 

A h y n = Y kl ~ a (Vn-k+i ~ 2y n - k + Vn-k-i) + n 1-a (Vi - Vo ), 

k=1 

because ai Q) = (n — 1 ) 1_ " — n l ~ a . □ 

Lemma 5. Suppose that y(x) is sufficiently differentiable function on [0 ,nh\ 

.. n—1 

—Ayn = h 2 ~ a Y kl ~ a y'n-k + (n/i) 1_a j/o.5 + O ( h 2 ) . 

fc=i 

Proof. From Lemma 4 


n—1 n—1 a 2 

^Ayn = X ^-“A^-fcW^Aiyo = h 2 ~ a Y k^-^Vn-k+n^h 

k =1 fc=l 


, 1— a i 1— a Ayo 


h 
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1 n —1 

—A h y„ = ft 2 -” V ft'-” (yl_ k + O (ft 2 )) + (nft)'-“ (j / 0 . 5 + O (ft 2 )) , 

k=l 

-j 72—1 / 72—1 

—A?/n = h 2 ~ a Y k 1 - a y"_k+(nh) 1 -°‘y , 0 5 +O ( h 2 ) (n/ft) 1- " + A" ^ fc 1_a 

fc=i V ?c=i 

The number ( nh) l ~ a ~ 0(1) is bounded. From d2J) we have 

72 — 1 

h 2 ~ a y kl ~ a ~ ^ 2_ “o (^ 2_Q ) ~ o(i). 

k=1 

Therefore 

.. n—1 

1 -Ayn = /A“ 51 k ^ a y'n—k + (n/i) 1_a yo. 5 + O (h 2 ) . 

k= 1 

□ 

Theorem 6. Lef y be a polynomial and x = nh. 

^AhV(x) = r (2 - a)y<«> + <(<*- l)j/"(i)ft 2 -“ + O (ft 2 ) . 

Proof. From Lemma 5 

1 re —1 

=ft 2 -“ £ e-X-* + (nft)'-Vo.5 + O (A 2 ) = 

fc=l 

72 

A 2 -” £ k 1 ~ a y"_ k - + U-X 5 + O (ft 2 ) , 

fc=l 

Then 

^Ai/(i) = (H ; - vD + O (^ 2 ) • 

From Claim 1 and Theorem 3 

= r ( 2 - a)J 2 ~ a y"(x) + + C(« - 1 )y'\x)h 2 ~ a + O ( h 2 ) , 

T(2 - a)r/ ( “)( x ) = r(2 - a)J 2 - a y"(x ) + ^(O)^ 1 - 0 . 
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Then 


£y" h = r(2 - a)y (a \x ) + C(« - l)y"(x)h 2 “ - x 1 a (y' 0 - + O (/r 2 ) , 

^A?/0) = T( 2 -a)?/ (Q) (a;)+C(a-l) 2 /' , (a;)/z 2 “"-x 1 " Q fy' 0 + ^ - y^+O (V) . 


By Taylor’s expansion 


y; + ^ = o (ft 2 ). 


Hence 


—A h y(x) = T(2 - a)y (a) (x) + C(« - 1 )y"(x)h 2 a + 0 ( h 2 ) . 


□ 

In Theorem 6 we determined the second order approximation for the 
Caputo derivative 


vL a> = 


T(2 — a)h a y T(2 — a) 
Corollary 7. Let y(x) be a polynomial. 


A h y n - ^ l \ &h 2 ~ a + O (h 2 ) . 


(7) 


y { ^ = 


T(2 - a)h° 


Y. 5 t ) Vn-k + 0(K 2 ), 


k =0 


where 5 j; Q) = afor 2 < k < n and 


^ = <?o a) ~ ((a - 1), 5{ a) = (t[ q) + 2C(a - 1), 5^ a) = a { 2 a) - ((a - 1). 

Proof. The second order backward difference approximation for the second 
derivative y" has accuracy 0(h). 


n Un 2y n _i + y n — 2 /o/m 

Vn = --+ °( h )- 


From approximation ([7]) 


!/£*> = 


-A h y, 


h 2 


C(« - 1) (Vn - 2y n _l + yn -2 


T(2 — a)h a * T(2 — a) 


h 2 


+ 0(h)) h 2 ~ a +0 (h 2 ) , 
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= --- -AhVn - _!L (Vn ~ 2 Vn-l + Vn- 2 ) + O ( K 2 ) , 


vr = 


r(2 - a)h° 


T(2 — a)h c 


<A = 


T(2 — a)h° 


^2 v^Vn-k - C(« - 1) (y n - 2y n _i + y n _ 2 ) j + O ( h 2 ) . 


K k =0 


□ 


Table 3: Error and order of approximation ([4j) for Caputo derivative of order 
a = 0.25 and y{x) = cos a; (left), y(x) = ln(a; + 1) (right) on [0,1]. 


h 

Error 

Order 

h 

Error 

Order 

0.05 

0.000081955 

2.30047 

0.05 

0.000029455 

2.31171 

0.025 

0.000017556 

2.22284 

0.025 

6.39 x 10~ 6 

2.20475 

0.0125 

3.95 x 10" 6 

2.15376 

0.0125 

1.46 x 10~ 6 

2.13162 

0.00625 

9.20 x 10~ 7 

2.10073 

0.00625 

3.44 x 10" 7 

2.08272 

0.003125 

2.20 x 10~ 7 

2.06368 

0.003125 

8.31 x 10~ 8 

2.05103 


Denote 

n 

B h y{x) = 22 s k 0) y( x - kh). 
k =0 

In Corollary 7 we showed that ((4)) is a second order approximation for the Ca¬ 
puto derivative of polynomials. Now we use the Weierstrass Approximation 
Theorem to extend the result to differentiable functions. 

Theorem 8. Let x = nh and y be a sufficiently differentiable function. 

y^ a \x) = ——-——— B h y(x) + O (h 2 ) . 

y V ; T(2 — a)h a J K J 

Proof. By the Weierstrass Approximation Theorem every continuous func¬ 
tion is a uniform limit of polynomials. Let e > 0 and p e (x) be a polynomial 
such that 

I y'(t) -p £ (t) | < e, 

for all t G [0, x\. Define 


Qe(t) = 2 /( 0 ) + 


Pe(f)df. 
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The function q e (t) is polynomial, and q'Jt) = p' e {t). We have that 


(y'(0 ~Pe(Q)d£ 


\y(t)-q e (t)\ = 
for all t G [0, x]. 

• r(l - a) | y (a \x) - < 7 e (a) (x)| 


< / \y'(Q-Pe(Q\d£< / ed£<xe, 


(x-Z) a 




< 


y^ a \x) — q^(x) < 


T(1 — a) J o 


(x- 0 " Q ^ = 


(x-O" 

ex 1_a 


dfl, 


T(2 — a) 


Therefore 


limg^fx) = y( a \x). 

e —>0 


Now we estimate Bh (y(x) — q e (x)). 


Y 4 “ } (Vn-k - Qe,n—k) 


| k=0 

We have that 


< Y l^l \ Vn ~ k ~ de,n-k\ < XeY I 5 , 


(«)l 

k 


k =0 


k =0 


E l 4 “’l < E ld“’l + 3|C(« -1)1 = 2- 3C(o - 1). 

k =0 k =0 


Hence 


I B h (y(x) - q e (x))\ < (2 - 3((a - l))xe, 

and 

lim B h q e {x) = B h y(x). 

e —>0 

• From Corollary 7 

= T(2-a)h° Bh<le ^ + ° ^ ' 

By letting e —>■ 0, we obtain 

y(a)(x) = r (2 - a )h? Bhv{x) + 0 (ft2 )' 

□ 
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Figure 1: Graphs of the Riemann zeta function on the intervals [—12, —1] 
and [—3, 0]. 


In Table 3 we compute the error and the numerical order of approximation 
(HD for the Caputo derivative of the functions y(x ) = cosx and y(x ) = 
ln(x + 1) on the interval [0,1], when a = 0.25. In Claim 9 and Lemma 10 we 
discuss the properties of the coefficients cr ^ and S 

Claim 9. Let 0 < a < 1 

- 0.1 < 4 a) < 0 . 

Proof. Denote 

a(a) = -4 1_a) = 2" +1 - 3" - 1. 

The function cr(a) has values cr(0) = u(l) = 0, and 


a'(a) = In 2.2" +1 - In 3.3". 


The first derivative of a (a) is zero when 


In 2.2" +1 = In .33", 



2. In 2 
In 3 


a = 


In 


2.1n2\ 
In 3 ) 


ln(3/2) 


« 0.5736. 


The function cr(a) is positive and has a maximum value <r(0.5736) ~ 0.0985 
on the interval [0,1]. □ 

The Riemann zeta function has zeroes at the negative even integers and 
is decreasing on the interval [—2,1], The value of ((a — 1) is negative, when 
a is between 0 and 1. Then Sq^ > 0 and < 0. From the properties of 
the coefficients of (]TJ) , the numbers S 1°° = aif 1 are negative, for n > 3. 
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Lemma 10. The number is positive when 0 < a < 1. 

Proof. From the definition of 5^ 

4 1_a) = 4 1_Q) - C(“«) = -o'(a) - C(-«) = z{a) - a(a). 

where z(a) = —£(—a). The function z(a) is decreasing on the interval [0,1] 
with values at the endpoints z(0) = 0.5 and z( 1) = 1/12 = 0.08333. The 
function cr(a) is increasing on the interval [0, 0.5736] and decreasing on the 
interval [0.5736,1]. 

Now we show that the minimum values of 5^ on the intervals [0, 0.8] 
and [0.8,1] are positive. 

min hi 1-0 '* > min z(a) — max cr(a), 
a€ [0,0.8] " ae [0,0.8] ae[0,0.8] 

min 4 1_a) > 2 ( 0 . 8 ) - a(0.5736) « 0.122 - 0.0985 = 0.0235. 

ae[0,0.8] 

and 

min 8^~ a ^ > min z(a) — max a (a ), 

<*€[0.8,1] <*€[0.8,1] <*€[0.8,1] 

min > z{\) - <7(0.8) ~ 0.083 - 0.074 = 0.009. 

<*€[0.8,1] 

Therefore the numbers 8^~°^ are positive when 0 < a < 1. □ 


5 Numerical Experiments 

In section 4 we showed that the approximation for the Caputo derivative 


< 7/(^0 


1 r(a) 

T(2 — a)h a k Vn ' 

has accuracy O (h 2 ) when n > 2. The numbers 8^ satisfy 

OO 

8 ( 0 a) > 0 , <5j a) < 0 , 8 ( 2 a) > 0 , 4 a) < 5^ < ■ ■ ■ < 4 a) < • • • < o, ^ 8[ o) = 0. 


k =0 
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In this section we compare the performance of the numerical solutions of the 
fractional relaxation and time-fractional subdiffusion equations using approx¬ 
imations ([T]) and (j3J) for Caputo derivative. From the Mean-Value theorem 
for the Caputo derivative 

y(h) - ?/(0) = h , y {a) (0), (0 <6 <h). 

1 (1 + a) 

The numbers T(1 + a) and T(2 — a) are between 0 and 1, when 0 < a < 1. 
Lemma 11. Let y be a sufficiently differentiable function on [0, h\ 

y(h) - 2 /( 0 ) - h a T(2 - a)y^\h) = O ( h 2 ) . (8) 


Proof. From the definition of the Caputo derivative 


y {a \h) = 


n J 


y'(0 


df. 


F(1 - a) J 0 (h-f) 

Expand the function yff) around f = 0 

2/(0 — y' (o) + fy" (o) + Off 2 ), 

rd-a), 


y' (o) 


Jo (h-O c 

We have that 


+ 2 /"( 0 ) 


£ 


(h~0 


;df + O ( h 2 ) 


f*h 


nh i _ h 1 -" r h 


:df = 


e 


df = 


I o (h-ty 
h 2 ~° 


-df. 


{ h ~ O a 1 - a ’ Jo { h ~ O a (1 — «)(2 — a ;) 


Then 


r (i - <«) = &(»)+ d-S-a ) ^”)+ 0 c-“)' 

F(2 - o)h a y (a \h ) = h (V(0) + ^2/"(0)) + O {h 3 ~ a ) , 

F(2 - a)h a y^(h) = h L'( 0) + ^"(0)) + O ( h 2 ) , 

F(2 - a)h a y^\h) = hy' Q) + O {h 2 ) = y(h ) - 2 /( 0 ) + O ( h 2 ) . 


□ 
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5.1 Numerical Solution of the Fractional Relaxation Equation 

The fractional relaxation equation (J2]) is an ordinary fractional differential 
equation with constant coefficients. The exact solution of the fractional re¬ 
laxation equation is determined with the Laplace transform method (lij . 
Numerical solutions of the relaxation equation are discussed in [19-21], In 
this section we compare the numerical solutions of the equation 

yW+y = F(t), (9) 


for approximations ([I]) and (dJ) of the Caputo derivative. When the solution 
y(t) of (jnj) is a continuously differentiable function, the initial condition y(0) 
is determined from the function F(t) by y(0) = F( 0). Let 


F(t) = 1 — At + 5t 2 


T(2 — a) 


t l ~ a + 


10 


, 2 -a 


r(3 


a 


Equation ([9]) has the solution 


y(t) = 1—4 1 + 5 1 2 , 


and the initial value y(0) = 1. Now we determine a second order numerical 
solution of ([9]) on the interval [0,1], using approximations ([4j) and (JHJ) for the 
Caputo derivative. 

Let h = 1/N, where N is a positive integer, and y n = y(x n ) = y(nh). In 
Lemma 11, we showed that (JBJ) 

V T{2-a)h° = yM( ' k ) + ° ( fe2 ~°) ' (10) 

Approximate the Caputo derivative y^ a \h) in equation ([9]) 

g^ + „(fc) = nfc)+o(ft^). 

Vi (1 + T(2 - a)h a ) =y 0 + T(2 - a)h a F l + O ( h 2 ) . 

Let {y fc }f =0 be an approximation for the exact solution y *. at the points 
Xk = kh. Set y 0 = y( 0) = 1. The value of yi is computed from the above 
approximation with accuracy O ( h 2 ) 

~ _Vo + r(2 - a)h a F 1 
Vl ~ 1 + T(2 - a)h a ' 
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The numbers y n , for n > 2, are computed from equation (JH]) by approximating 
the Caputo derivative y^ with (J3J • 

1 n 

T(2 - a )h a ^ ^ )?/n_fc + Vn = Fn + 0 (^ 2 ) ’ 

n 

y n + T(2 - a)fr“) = T(2 - a)h a F n - £ + O (/r 2+ “) . 

k =1 

The numerical solution {^fc}^ =0 , for 2 < n < N, is computed explicitly with 

= C + T(2 - a)h° ( F<2 “ aWF " ~ ' (U) 

Similarly, we obtain an explicit formula for the numerical solution \ykj k _ Q 
of equation ([9]), by approximating the Caputo derivative yn with (JTJ) 

*• = i + rp-W^ ( F(2 “ a,haK ~ t ' <12) 

Numerical solution (HU) converges faster to the solution of the fractional 
relaxation equation, because it has a second order accuracy O ( h 2 ), and the 
accuracy of numerical solution (1121) is O {h 2 ~ a ). 

Figure 2: Graph of the exact solution of equation (JU]) and numerical solutions 
(fTTlb black. and (TT2l) -red. for h = 0.1 and a = 0.8. 
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Table 4: Maximum error and order of numerical solutions (fl2l) and (ITT]) for 
equation ([9]) on the interval [0,1], when a = 0.8. 


h 

Error 

Order 

0.05 

0.0628014 

1.17381 

0.025 

0.0275997 

1.19262 

0.0125 

0.0120751 

1.19262 

0.00625 

0.0052704 

1.19603 

0.003125 

0.0022975 

1.19785 


h 

Error 

Order 

0.05 

0.0081544 

1.85708 

0.025 

0.0021629 

1.91461 

0.0125 

0.0005599 

1.94979 

0.00625 

0.0001428 

1.97076 

0.003125 

0.0000361 

1.98307 


5.2 Numerical Solution of the Fractional Subdiffusion Equation 

The time-fractional fractional subdiffusion equation is obtained from the heat 
transfer equation by replacing the time derivative with a fractional derivative 
of order a, where 0 < a < 1 

d a u(x,t) _ d 2 u(x,t ) 
dt a dx 2 

with initial and boundary conditions 



u(x, 0) = u 0 (x), u(Q,t) = u L (t), u(l,t) = u R (t). 


Numerical solutions of the fractional subdiffusion equation are discussed in 


19), 21, 24, |34]. In this section we determine the numerical solutions (fT5l) 


and (1TBD for the fractional subdiffusion equation obtained by approximating 
the Caputo derivative with (P and (0]) on the region 

(x,t) e [0,1] x [0,1]. 

Let h — 1/N, t — 1/M, where M and N are positive integers, and Q be 
a grid on the square [0,1] x [0,1] 

Q = {( nh , mr)||l < n < N, 1 < m < M} . 

Denote by u™ and F™ the values of the functions u(x,t ) and F(x,t ) on Q 

u™ = u(nh, mr ), F™ = F(nh, mr). 


By approximating the values of the Caputo derivative in the time direction 
at the points (nh, r) using (TTOT) and using a central difference approximation 
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for the second derivative in the space direction we obtain 


ui - w° u*_i - 2id, + u\ +x /19 2 -«\ 

r°f(2 - a) = // -^±i + F( n fc, t) + O (/r + r 2 ") . 

Let 

, = r( 2 -a)L. 

The solution of the fractional subdiffusion equation satisfies 

~wl-i + (1 + 2r?X - r?< +1 = + T(2 - a)r a F(7ih, t) + O (r a h 2 + r 2 ) . 


Let 17™ be the numerical solution of the fractional subdiffusion equation on 
the grid Q. The numbers 17™ are approximations for the values of the solution 
«™ = u(nh,mT). The numbers 17° are computed from the initial condition 
17° = u 0 (nh). The numbers 17* are approximations for the solution of (1151) 
at time t — r. We compute the numbers 17* implicitly from the equations 

-riULi + (1 + 2r,)t/‘ - > lK +l =U° n + T(2 - a) T ‘Fl 
where the values of Uq and U ^ are determined from the boundary conditions 

Uq = u l {t ), 17* = u r (t). 


The numbers 17* are computed with the linear system (k — 2, ■ ■ • , N — 2) 

( (1 + 2,rj)Ul - rjU^ = u 0 (h) + rju L (r) + T(2 - a)r a Fl 
{ ~vUl_i + (1 + 2rj)Ul - r)U£ +1 = u Q (kh) + T(2 - a)r Q F fc * 

{ -vUh -2 + (! + 2 v)Uh-i = Uq((N ~ 1 )h) + VUr(t) + T(2 - a)r a F^_ v 


Let K be a tridiagonal matrix of dimension 77 — 1 with values 1 + 2 rj on the 
main diagonal, and —rj on the diagonals above and below the main diagonal. 


and 17™ 
system 


/l + 2 1 ] —r) 0 0 0 \ 

—i/ 1 + 2ij —rj 0 0 

I< 5 = 0 -i) l + 2i| -i) 0 

0 0 —rj 1 + 2 rj —rj 

^ 0 0 0 — 7] 1 + 2rjJ 

(17™, 17™, • • • ,17™_ 1 ). The vector 17 1 is solution of the linear 
KU 1 = R 1 + 7]R 2 , (14) 
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where R\ and R 2 are the column vectors 

Ri= [u 0 (kh) + T(2-a)r a F^, 

R2 = [■ u l {t ), 0 , • • • , 0 , u r {t)] t . 

We determined a second order approximation U 1 for the solution of the 
fractional subdiffusion equation, on the first layer of Q, as a solution of the 
linear system (I14p . When m > 2 we discretize the Caputo derivative with 
equation (115]) with the second order approximation ([1]) 


1 


t°T(2 — a) 


\ ^ r ( Q ) m —k U rt— 1 — . . ( 2 2\ 

—— 2^ u n = - jj -— + F(nh, mr) + O ( h 2 + r 2 ) . 

a ' k=0 1 


The values of the numerical solution U™ are determined from the equations 


-nUT-1 + (<So° + 2ri)ic - ’Kl 1 = -J2 s P u ”~ t + r ( 2 - “K-C, 

fe=l 


and the boundary conditions 

U™ = u L (mT), U™ = u R (mr). 
The vector U m is a solution of the linear system 


(K — ((a — l)I)U m — Ri + r]R 2 , 


(15) 


where R\ and R 2 are the column vectors 


Ri = 


-£4 a| tC'‘ + r( 2 

. k= 1 


a r 


N-l 


n =1 


R 2 = [u L (mr), 0, • • • ,0, u R (mT)) T . 

The numerical solution j[/ 2 , • • • , U M }, using approximation (J4]) for the Ca¬ 
puto derivative, is computed with linear systems (fT5lh Similarly we deter¬ 
mine the numerical solution {V 2 , ■ ■ ■ , V M } for approximation ([TJ) with linear 
system (1T6l) and hrst layer V 1 — U 1 . 

The numerical solution V m is computed with the linear system 


KV m = R 1 + V R 2 


(16) 
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where R\ and i? 2 are the vectors 


Ri = 


2^ a k 

k=\ 


>v r 


+ r (2 - a)r“F; 


N-l 


-I n= 1 


R'2 = [u L (mT), 0, • • • ,0, u R (mT)) T . 

Numerical solution ([1311 has accuracy O ( h 2 ) and the accuracy of (HUD is 
O (h 2 ~ a ). When 


lOt 


2 -a 


41 


r(3 - a) r(2 -a))' 


F(x, t ) = 2(1 — 3a;)(5 1 2 — 4t + 1) + a: 2 (l — x ) 

the fractional sub-diffusion equation 

d a u(x,t) d 2 u(x,t ) . . . . r n r n 

at* -—^r L + FI ' x ' t ^ ^ *) 6 l°-1] x [o, i] 

with initial and boundary conditions 

u(x, 0) = x 2 (l — x), u(0 , t) = u( 1, t) = 0, 


1 —a 


(17) 


has solution 

u(x, t ) = ic 2 (l — x)(l — 4t + 5t 2 ). 

The maximal error and numerical order of numerical solutions (TT6|) and (fl5l) 
for t = h and r = h /2 at time t = 1 for the fractional subdiffusion equation 
(TT7T) are given in Table 5 and Table 6. 


Table 5: Maximum error and order of numerical solutions (fl6l) and (fT5]l for 
equation (1T7T) when a = 0.6 and t — h, at time t — 1. 


h (r — h ) 

Error 

Order 

0.05 

0.00051794 

1.37686 

0.025 

0.00019766 

1.38974 

0.0125 

0.00007530 

1.39222 

0.00625 

0.00002864 

1.39467 

0.003125 

0.00001087 

1.39657 


r (r = h ) 

Error 

Order 

0.05 

0.00001170 

1.93892 

0.025 

2.99 x 10" 6 

1.96593 

0.0125 

7.62 x 10~ 7 

1.97559 

0.00625 

1.93 x 10" 7 

1.98175 

0.003125 

4.87 x 10~ 8 

1.98648 
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Table 6: Maximum error and order of numerical solutions (1161) and (1T5]) for 
equation (ITTT) when a = 0.4 and r = 0.5 h, at time t — 1 . 


h (h = 2 t ) 

Error 

Order 

0.05 

0.00006172 

1.56262 

0.025 

0.00002069 

1.57697 

0.0125 

6.91 x 10" 6 

1.58139 

0.00625 

2.30 x 10" 6 

1.58597 

0.003125 

7.65 x 10" 7 

1.58946 


h (h = 2 t) 

Error 

Order 

0.025 

4.41 x 10" 6 

1.98507 

0.0125 

1.11 x 10" 6 

1.99521 

0.00625 

2.78 x 10~ 7 

1.99589 

0.00625 

6.95 x 10" 8 

1.99701 

0.003125 

1.74 x 10" 8 

1.99811 


In numerical solutions (fl~5li and (TT6|) . we use (fTOl) to obtain a second order 
approximation for the solution of the fractional subdiffusion equation on the 
first layer of Q. Another way to determine a second order approximation for 
the solution at time t — r is to compute the partial derivative u t (x,t ) at 
time t — 0 and approximate the solution u(x,t) with a second order Taylor 
expansion. The function u(x,t ) satisfies 


d a u(x,t ) d 2 u(x,t ) 

dt a dx 2 

Apply fractional differentiation of order 1 — 


+ F(x,t). 


d 1 “ d a u(x,t) d 1 “ d 2 u(x,t) d 1 a F(x,t) 

dt l ~ a dt a dt 1 ~ a dx 2 <9t 1-Q 

9 1_Q d 2 u(x,t) d 1 ~ a F( y x,t) 

M x ^) = dfl . a dx 2 + dfl „ a • 


Set t — 0 


u t (x, 0) = 


d 1 " d 2 u(x,t) 


dt 1 - 


dx 2 


+ 


d 1 a F(x,t) 


t =o 


dt 


1 —a 


t =0 


When the solution u(x,t) is a sufficiently smooth function 


we obtain 


d 1 a d 2 u(x,t ) 
dt 1-Q dx 2 


Ut(x, 0) = 


d 1 a F(x,t) 


dt 1 ~ c 


t =o 
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The values of the solution at time t — r are approximated using the second 
order Taylor expansion 


u(x, t ) = u(x, 0) + ru t (x, 0) + O (r 2 ) 
In the fractional subdiffusion equation (H71i 

lot 2 -" 


F(x, t ) = 2(1 — 3x)(5 t 2 — At + 1) + x 2 (l — x) 

10t 1+ “ 4 t a 


At 


1 —a 


- 2(1 - 3 .) 


T(2 + a) T(l + a) 


T(3 - a) T(2 -a))' 
+ x 2 (l — x) (lOt — 4). 


The partial derivative Ut(x,t) at time t — 0 has values 

d 1 ~ a F(x, t ) 


u t (x, 0) = 


dF- c 


= —4x 2 (l — x). 


t =o 


Then 


u(x, r ) = u(x, 0) + ru t (x, 0) + O (r 2 ) 

u(x, t ) = x 2 (l — x) — 4ra; 2 (l — x) — x 2 (l — x)(l — At) + O (r 2 ) 

We obtain the second order approximation for the solution of equation (fT71) 
on the first layer on the grid Q. 

U,l = {nh) 2 {l - n/0(l - 4r), (n = 1, 2, • • • , N - 1) (18) 


Table 7: Maximum error and order of numerical solutions (I15p and (fIBli with 
approximation (1T8|) for the solution of equation (TT7l) on the first layer of Q , 
at time t — 1 when a = 0.6 and r = h. 


h (r = h) 

Error 

Order 

0.05 

0.00051282 

1.35060 

0.025 

0.00019690 

1.38100 

0.0125 

0.00007518 

1.38896 

0.00625 

0.00002862 

1.39337 

0.003125 

0.00001088 

1.39601 


r (r = h) 

Error 

Order 

0.05 

0.00001730 

2.28453 

0.025 

3.85 x 10" 6 

2.16657 

0.0125 

9.02 x 10~ 7 

2.09523 

0.00625 

2.17 x 10~ 7 

2.05667 

0.003125 

5.29 x 10~ 8 

2.03510 
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6 Conclusion 


In section 4 we compared the numerical solutions of the ordinary fractional 
relaxation equation and the partial fractional subdiffusion equation using 
approximations ([I]) and (j4j) for the Caputo derivative. The higher accuracy of 
approximation (J3J) results in a noticeable improvement in the performance of 
the numerical solutions. Numerical experiments suggest that the numerical 
solutions converge to the exact solutions of the fractional relaxation and 
subdiffusion equations for all a between 0 and 1. We are going to work on a 
proof for the convergence of the numerical solutions discussed in section 4. 
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